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Preface 

This module is directed toward intermediate undergraduate students 
of the ecological sciences who are familiar with differential calculus and 
algebraic techniques. No skill in computer or calculator programming is 
required. Some of the mathematical functions used are fabricated to illustrate 
a mathematical technique or concept. However, each function's behavior is always 
comparable to the characteristics of the ecosystem it represents. Many topics 
are introduced and extended in the problem sets. Thus the problems and 
computer exercises do not merely repeat the text but form an important 
supplement. 
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INTRODUCTION 

Most applications of calculus to ecological problems involve the 
determination of specific relationships between measured quantities. 
Frequently, the obvious relations involve rates of change of the quantities 
of interest, and not the quantities themselves. In mathematical terms, we 
often can determine equations involving the derivatives of functions instead 
of the functions themselves* We then need to use an inverse operation on the 
derivative to determine the desired function. This operation is called 
integration and its field of mathematical study is called integral calculus. 

For example, we may wish to know the blood level of a toxic material 
that is absorbed through the skin during a certain time period. However, the 
mathematical model may be based on the rate of absorption through the skin and 
rate of excretion via the urine. Thus we wish to know the magnitude of a 
quantity (body level of pollutant) when we know only the rate of change of 
that quantity (rate of increase by absorption and rate of decrease by excretion). 



FUNDAMENTALS 

The inverse operation to differentiation is called antidif ferentiation 
or, more coimnonly, integration . In models of exponential growth we assume 
that the growth rate is proportional to the population size. 

f-kK. (1) 
where N is the population size, t is time ^nd k is a proportionality constant. 
We could solve for N as a function of time by integrating the derivative. 
However, we know that the exponential function is the only function which 



2 



equals its derivative. 



dt 



Thus, we can modify the function to give the solution 

kt 



N = N e 
o 

kt 

since differentiating N^e gives eqn. 1. 

d(N e^*') . 

^ = 2 kN e'''' - kN 

dt dt o 

If we write eqn. (1) as 

^ «= kN e'^^ (2) 
dt o 

then the "antidif f erentiation" becomes more obvious: find N so that its 

kt 

derivative equals kN e . If we define 

o 

g(t) « ^NqS^^* f (t) = N^e^^ 
then eqn. (2) becomes 

g(t) = df/dt. 

Thus g(t) is the derivative of f(t) and, conversely, f(t) is an antiderivative 
of g(t). But there is a slight problem. The antiderivative is not unique. 
If we define 

x(t) « N e^^ + 10, y(t) « N e^^ + 354 
o o 

then x(t) and y(t) are also antiderivatives of g(t), as can be checked by 
differentiating: dx/dt « g(t), dy/dt = g(t). One interpretation of this 
is that the graphs of x(t) and y(t) have the same slope (derivative) for 
any given value of t. Since the functions differ by only a constant, we can 
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write the general form of the antiderivative as f (t) + C where C can be any 
constant. We usually write the antiderivative as the indefinite Integral 

g(t)dt « f(t) + C (3) 

and call C the integration constant . Recall that the term "dt" identifies 
the variable of integration just as it does the variable of differentiation 
in the derivative df/dt. 



THE DEFINITE INTEGRAL 



The integration constant does not appear if the integral is a definite 
integral , written as 



A = 



g(x)dx 



where a,g are called the limits of the integral* In this case, the integral is 
determined by evaluating the antiderivative of g(x) at the values x « b, x « a and 
subtracting. Thus, as with eqn. (3), if 

df f V 
d5F= 



then 



g(x)dx = [f (x)] 



= f(b) - f(a), 



Note that the integration constant is not written since it is eliminated 
through subtraction. If g(x) = Zx^, then |g(x)dx - f (x) + C - 2x^/3 + C and 
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g(x)dx - 2x^dx - (2bV3 + C) - (2aV3 + C) 

- 2bV3 - 2aV3 + C - C 
= f(b) - f(a). 



Area Under the Curve 

Whereas the derivative can be used to represent the slope of a curve, the 
definite integral can represent the area under the curve (specifically, the 
area between the curve and the horizontal axis). For a curve given by y - g(x), 
the area (A) under the curve from x -= a to x = b is written (see fig. 1). 



g(x) dx 




a b 
Fig 1. Area under the curve g(x). 

This use of the definite integral has strong intuitive appeal. Consider the 
case where g(x) = 1 as in figure 2. 
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a b 
Fig. Area of rectangle. 

The area under g(x) is then 



A = 



g(x)dx 



2 dx « 2 



dx •= 2(b-a). 



'a 



Thus 



dx represents the width and g(x) is the height. Thus the definite 



integral as area has intuitive meaning: area « height x width. In general. 



the same visual identification applies: 

rb 



g(x) dx 



width ^ 
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Example 1 

Diffusion of a gas across a membrane is often described by Pick's Law which 
states that the rate of transport is proportional to the product of the surface 
area of the membrane and the concentration gradient. Thus knowledge of the surface 
area is Important. Certain leaves have nearly parabolic edges (fig. 3). The 
area of the top surface available for transport of water is then found by integrating 

10 
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the parabola function. If the dimensions of the leaf are length * 2a, width « 2b, 
then the parabola describing one edge is written 

y(x) « -bxVa^ + b. 




Fig. 3. Leaf with parabolic edge contours. 



The area of half the leaf is 

(^bx^/a^ + b)dx 



A 
2 



-a 



i.e., the area between the x-axis and the upper curve. 



Thus we calculate 



A = 2 



(-bxVa^ + b)dx 



-a 



« 2[-(b/a2)xV3 + bx]^ 

■= 2[(-ab/3 + ab) - (ab/3 - ab)] 

= 8ab/3. 
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Improper Inte^^rals 

Occasionally we are interested in long term behavior of a measured quantity. 
For example, when raw sewage and other chemical pollutants are continually dumped 
into a lake, one effect is the rapid increase in the number of microorganisms. 
The resultant high level of organic oxidation from metabolism of the sewage can 
lead to extreme oxygen depletion of the lake, with obvious detrimental effects 
on other aquatic life (Dugan 1972). If we know something about the rate of oxygen 
consumption by the microorganisms, then we obtain the amount of oxygen consumed 
during a time period T by integrating the rate over that time period, i.e., using 
a definite integral with limits 0,T. We then estimate the maximum oxygen 
depletion by integrating over an infinite time interval. Such a definite integral 
with at least one infinite limit is called an improper integral . 
Example 2 

Assume that self-inhibition by the microorganism population causes the 
oxygen depletion to taper off as time becomes large. One model might be 

^ -t^. -t 

— - te + e (4) 

where C represents the quantity of oxygen consumed. A graph of dC/dt looks like 
figure 4. 
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Fig. A Rate of oxygen consumption. 



The total amount of oxygen consumed by time T is then found by integrating eqn. (4), 

fT 



C(T) 



(te~*^ + e~*^)dt 



The maximum depletion is then found by letting T increase to 

fT 



C(oo) « lim 

X-XX) 



(te~*^ + e~*^)dt 



With most nicely behaved models, we can evaluate C(co) directly. The above integral 
is then written 



C(«.) 



(te~*^ + e'"*^)dt 



(5) 
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and Is evaluated as with previous definite integrals. 

Some functions are not "nicely behavi^:''" and certain integrals cannot be 
evaluated directly. If we desired to evaluate the circumference C of a tube with 
circular cross-section by using an integral (instead of the well<-knovn formula 
C » 2TTr) we can encounter difficulties. For one-quarter of the circle (fig. 6) 
the arc-length is given by (Schwartz 197 A, p. 63A) 



rr 



L = 



dx 



where r is the radius. The integrand is infinite at x "r (fig. 7). However, we 
can evaluate the integral by taking the limit of another integral: 



lim 



dx 



-1 

= lim [r sin (x/r)] 
a-+T 0 

= lim [r sin"-'-(a/r)] - tti/I. 
a->-r 
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METHODS OF INTEGRATION 



Since there are no foolproof formulae that can be used for all integrals 
(as there are for derivatives), the following integration methods all attempt 
to change a difficult integral into a simpler one. All examples use definite 
integrals. 

Substitution 

This method substitutes each part of a definite integral with a counterpart 

so that the result is kept the same. With the substitution u « h(x), the integral 

'b 

f (x)dx 

becomes transformed into 

'h(b) 
g(u)du 

h(a) 

where g(u)du = f (x)dx. Thus the integrand, dx, and the limits have been transformed 
into equivalent counterparts. We obtain g(u) by finding the inverse function 
h ^(u) so that 

X = h"-'-(u) 

d[h^^(u)] . 
dx « — — , ^ du 
du 

g(u) du = [f (x)][dx] 

g(u) = [f(h-\u))] [^^^^]. 
Often, h(x) appears explicitly in the integrand so that u is substituted directly. 
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Example 3 

The problem presented above on the oxygen depletion due to sewage dumping 
can now be solved. The integral in eqn. (5) is first separated. 



(te"*^ + e"^)dt - 



te dt + 



|e"^dt 



(6) 



The antiderivative of te- is not obvious, so we substitute a new function into 
the integral. Let u «= t^. Then the differential du is 
du = (du/dt)dt « 2tdt 

so that 

du 

The limits remain, since u « 0 when t = 0 and u « " when t = ~. Direct 
substitution then gives 



te dt 
0 



e"^ (tdt) 



0 

mOO 



-u >duv 1 
e (— ) = 2 



e du 



0 



l[.e-]" =i[o- (-1)1 
^ 0 



3, 
2 



Since the second term in eqn. (6) is now obvious ( 
we obtain the solution to eqn. (5) of 

CC") = 1/2 + 1 = 3/2. 



e du 



e~*^dt = 1), 
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Integration by Parts 



This method utilizes a relation from differential calculus concerning the 
total differential. Recall that the differential of a product of two functions 
u(x), v(x) can be written 

d(uv) = udv + vdu. 



If we evaluate ant Iderlvat Ives » we obtain 



uv 



udv + 



vdu. 



Rearrangement of eqn. (7) gives the Integration by parts formula 



udv = uv - 



vdu 



(7) 



With definite integrals, we usually write this formula as 

b b 

b 



[u(x) ^]dx = [uv] - 
a 



[v(x) f ]dx 



a 



(8) 



Again, the goal Is to change a difficult Integral, the left side of eqn. (8), 
Into a simpler one, the right side of eqn. (8). An example Is presented later. 



Partial Fractions 



When the Integrand Is a ratio of two polynomials. It often can be decomposed 
Into a sum of simpler terms. Only the case of nonrepeated linear terms In the 
denominator Is treated here. For more complicated cases In an ecological setting, 
see Clow and Urquhart (1974) p. 559. 

18 
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The integrand is assumed to be of the form F(x)/G(x) where F(x) and G(x) 
are polynomial functions and where G(x) is the product of linear factors. For 
example, 

G(x) « (1 + 2x)(2 + 2x)(l + x) 

is a polynomial composed of factors linear in x. The partial fraction technique 
replaces the single rational expression by a sum of terms where each denomination 
is one of the linear factors of G(x), If we have two linear factors, 

G(x) = A(x)B(x) 

then a partial fractions decomposition gives 

F(x) _ F(x) ^1 , ^2 

G(x) ^ A(x)B(x) " A(x) B(x)' 



where and are constants. Multiplication by A(x)B(x) gives 



F(x) = C^B(x) + C2A(x) 



(9) 



Let r^,r2 be zeros of A(x), B(x), respectively, i.e. A(r^) = B(r2) = 0. Then, 
with X = r, eqn. (0) is 



F(r^) = B(r^) 



and with x « r^^ eqn. (9) becomes 



F(r2) = C2 A(r2) 



so that ^^9^2 easily determined, 



Example 4 

The logistic growth model for animal populations is represented as a 
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differential equation 



dN 



^ » rN(l - N/K), r,K « constant. 

Separating variables (see the section, Differential Equations) gives, using 
the differentials dN and dt. 



dN 



N(l - N/K) 



rdt 



which is integrated to yield the equation 

dN 



N(l - N/K) 



rdt 



(10) 



The right side of eqn. (10) is easily integrated. The left side, however, must 
be reworked. First rewrite the integral as 



N(l - N/K) N(K - N) 



Now expand in partial fractions as 



K ^1 ^ ^2 



N(K - N) N K-N 



Mulitplying both sides by N(K - N) gives 



K = a^(K - N) + a^N. 



(11) 



Since (11) must hold for all values of N (Clow and Urquhart 1974, p. 562), then 
setting N = 0 gives 



and N " K gives 



so that 



a^^ = 1 



a^^ 1 



K 



4+ 1 



N(K - N) N K-N * 
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Then eqn. (10) is 



dN 


» 

dN 


N 


K-N 



which integrates to 

iln N - iln (K-N) « rt + C. 

APPLICATIONS OF DEFINITE INTEGRALS 

Direct applications of integrals generally fall into discrete categories 
in contrast to applications of derivatives which usually are based on slopes. 
The first group discussed below uses the integral as the accumulation of changes 
in the function. The second category uses the integral as an area or generalized 
volume. The last application is more mathematical, although it actually relies 
on the accumulation concept, and uses the integral to estimate the error in a 
given approximation. 

Accumulation of Changes in the Function 

The integral as a total accumulation has been presented before in example 2 
on oxygen depletion. This use of the integral is actually fairly intuitive. Let 
us call our quantity of interest F(x). Then F' (x) » dF/dx is certainly the rate 
of change of F(x) and F(x) is certainly the antiderivative of F'(x). Then 
integrating the rate of change of F gives the total change in F. 

b 

F'(x)dx = F(x)] « F(b) - F(a) . 
a 

^a 

Thus the definite integral 
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b 

F'(x)dx 

Is the total change In F(x) as x changes from a to 
Average Change 

The average change in F(x) is then found by dividing by the change in x, 
since the average is the change in F per unit change in x. Note that this 
formula can be shown graphically as the average height of the function. For 
a given curve, the area under the curve equals the average height multiplied 
by the width. Thus the average height y of a curve y « f(x) is the area A 
divided by the width. 

b 

f (x)dx » (b-a)y 

a 



- „ ^ 1 
y " b-a .b-a 



b 

f (x)dx 

a 



Example 5 

An interesting example is a study (Fisher 1963 in Warren 1971, pp. 161-163) 
of the effects of dissolved oxygen content and food ration on the growth rate of 
Coho salmon. The data appear in figure 7. The upper curve is well approximated by 

y - 7.3(x + 3.5)e"-°^^ (12) 

The lower curve is the straight line 

y « 28 

where y » growth rate, x « dissolved oxygen. A simple comparison of the effect 
of diet (restricted vs. unrestricted ration) on growth rate is to compare the 
average growth rates (y) for the two diets. Since the lower curve has a constant 
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60- 



>> 

^ 50- 

to 



to 
a 



J-i 3 0- 



S 20- 



unrestricted ration 




(3.0-18.0 tng/1) 
(3.0-9.5 tng/1) 

restricted ration 



— □ — □ — cr 



1 0- 



10 



Dissolved oxygen (mg/l) 



20 



Figure 7. Relationships between dissolved oxygen concentration and growth 
rate of juvenile coho salmon when food was unlimited and when it 
was limited. Arrows indicate growth of fish when held at oxygen 
concentrations fluctuating diurnally between levels specified. 
Data of Fisher (1963), in Warren, 1971, p. 162. 
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growth rate, we must have y - 28. 
- 1 



27 



30 



28dx - ■^[28x]^ - -^[(28) (30) - (28) (3)] 



- jj (28) (27) - 28 



For the upper curve. 



- 1 



y = 



27 



30 



7.3(x + 3.5)e"°'^^dx 



7.3 

27 



= .95 



30 30 

- _ -.0.5x, , 7.3 -.05x^ 
3.5e dx + ~27~ ^'^ 



30 30 

-.05x, , -.05x. 
e dx + .27 xe dx 



-.05x, 1 -.05x. 

e dx = .95[- e ] 



30 



The first integral is the same form as in previous examples. 

|-30 

.95 

= 'l9(-e--°^(30) ^^-.05(3)) ^ ^2.07 

Rather than finding the second integral in tables, we will evaluate it using 
integration by parts . 
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With dv(x) representing the differential of v(x), the second integral can 
be written 



f30 



30 



-.05x^ 

xe dx 



/ ^ J / \ -.05x, 

vhere u(x) » x, dv(x) » e dx 



u(x) dv(x) 



Recall that integration by parts uses the formula 
judv » uv - |vdu 



Since du = dx, and v «= 
f30 

-.05x 

xe 
3 



dv = -e"-°^^/.05, we get 



dx «= [ 



-xe 
.05 



-.05x 30 
] 



30 



-.05x 



.05 



dx 



1 

= .05 



- . -.05x 30 

[-30e~^-^ + 3e~-"] - ] 



, . -.15 -1.5 

20[-30e~^-^ + 3e~-^^ + ^-^^ - -—r" 1 



.05 



.05 



= 173.04 



Thus 

7 = 12.07 + .27(173.04) = 58.8. 
In sximmary, we have the averages: 

ration average growth rate 

restricted 29.3 

unrestricted 58.8 
It is somewhat surprising that the average unrestricted ration rate is over twice 
that of the rate for the restricted ration. The data is deceptive visually due 
to the close values near x = 3 mg/.H and the distorted logarithmic scale. 
Distance 

Velocity is defined as the rate of change of position. Since the distance 
covered is the total change in position, it must equal the integral of the velocii 
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For a free-falling object, the velocity is given by v(t) ■ g*t where g is 

the acceleration due to gravity and t is time elapsed. The distance S covered 

1 , 

after T seconds is given by S - •5- gT which is merely the integral of velocity. 



v(t)dt 



(gt)dt 



t dt - g[tV2] 



Volumes 

Volumes and areas of complicated regions are also evaluated using the 
definite Integral. Previously, the area under a curve was bounded by three 
straight perpendicular lines. When the bottom is not the base axis, the 
integration is still simple. For a region shown in figure 8 the area is the 
difference between the area under curve f and the area under curve g . Thus 



f (x)dx - 



g(x)dx 



a 

b 



a 



[f(x)-g(x)]dx 



Note that [f (x)-g(x)] is merely the height of the region at the point x, so 
the height times width interpretation is still applicable. 



26 



22 




f(x) 



g(x) 



X 



a 



b 



Fig. 8. Area between two curves. 



Surface Area of Revolution 

When the region Is not planar, the evaluation of Its area must take Into 
account the changes In the third dimension. If the surface is obtained by revolvlr 
a curve around a straight line, the evaluation needs only a single Integral. The 
following example Illustrates the method. ^ 
Example 6 

One study of temperature regulation in mammals requires knowledge of the 
surface area exposed to the sun. The model views the torso of the animal as 
symmetric with respect to a longitudinal axis. Each vertical cross-section is 
then a circle. The simplest such approximation is a cylinder: 



The cylinder can be described by revolving a straight line around the axis, as 
in figure 9 • 




torso 




tail 
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Unrolling the surface gives a rectangle wh ose width equals jt he^-clr-cumference-pf ^ 
the circular end face of the cylinder. This surface area of revolution equals 
the line length multiplied by the width, thus (fig. 10) 

A « l-nr^i 




Fig. 10. Surface of unrolled cylinder. 



As with the area under a curve, the general formula for a surface area of 
revolution must be intuitive, i.e., must visually appear as length times circum- 
ference. Let the torso have a profile of varying radius: 

10 



28 



2A 





Ir(x) ^ 











Assume the line to be revolved is graphed as follows, 




and is represented by r « a - bx^, a and b positive constants. The radius then 
changes with x, and the integral must be used: 



= I 2TTr(x)d3 



I circumfer 



circumference 



In this example, let x- ■ .5m, a « .28, b - .2A. The total surface area is 



.5 



27r(.28 - .2Ax2)dx 



-.5 



- 27r[.28x - x'] 



.5 



-.5 



- 1.63 m^ 

Note that any function will work in the formula, as long as the area desired is 
a surface area of revolution. The only problem might be in using a function which 
Is difficult to Integrate. 
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Volume of Revolution 



The volume of a solid generated by revolving a curve around an axis can 
be derived as an intuitive extension of the surface area of revolution. The 
area and volume formulae for the cylinder and the general revolved solid 
(figure 11) are seen to be analogous. 



Surface area 



Volume 



Cylinder 



27rr^£ 



Trr^il 



Revolved solid 



27rf(y)dy 



7r[f(y)]2 dy 





Fig. 11. Solids of revolution. 



Now we develop the general volume formula by expressing the integral as a 
limit of sums of pieces of the solid. Consider a curve z ■ f(y) and the region R 
under the curve (figure 12). We revolve the region R around the y-axis (figure 13 
to obtain the solid. 
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b y 




Fig. 12. The function f(y). 



Fig, 13. The solid obtained by 
revolving f(y). 



Suppose we divide the interval [a,b] into many subintervals, each of width 
dy. Then, if dy is sufficiently small, the area of the subregion is well 
approximated by a rectangle of width dy and height f(y^), as figure 14 indicates, 
By revolving about the y-axis, we sweep out a circular slab with radius f (y^) 
and thickness dy (figure 15). 





Fig. 14. The area increment dy. 



Fig. 15. The volume element obtained 
by revolving dy. 
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The volume of each slab is it (radius)^ (length) « 7Tf(y^)^dy. Thus we have volume V 
of the solid as 



V « lim ^ l TTf (y )2 dy, 
dy-K) i»l 



TTf (y)^dy. 



N - (b-a)/dy 



General Surface Areas 

When the surface is more irregular and is not axially symmetric, its area 
can still be found. The surface must now, however, be described by a three- 
dimensional function which gives the height as a function of the length and 
width coordinate: z « f(x,y). 




The dependence on two variables requires two integrals, and the method used is 
called double integration ^ 

Given a function of two variables, say z " f (x,y), we can write a double 
integral of z over a region R as: 



f (x,y)dA 



R 
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with A representing the area coordinates from the region R. This integral is 
more often evaluated by writing it as an iterated double integral ; 



fb 



a 



h(x) 



g(x) 



f(x,y)dy dx 



h(x) 



f (x,y)dy 



dx 



Evaluation of the "inner" integral yields a function of x, which becomes the 
integrand of the "outer" integral.* 

When the surface is described by z « f(x,y), its area is found using an 
iterated integral. The limits of integration are found by projecting the boundary 
of the surface onto the x,y plane. The formula for the surface area is** 



A = 



h(x) 



g(x) 



and is best illustrated by example. 
Example 7 

If the animal is again considered-to look like a cylinder, one improvement 
would be to account for the neck rising at an angle from the shoulder. To keep 
the calculations simple, we assume the neck rises vertically, and is also 
cylindrical (figure 16). 




Fig. 16. Half of the upper surface. 



♦Note that the order of integration can be reversed when f (x,y) is cont.t:aous in x an 
**See Ayres (196A), p. 319. 
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We determine the area of the torso by subtracting the "back" area inside the 
vertical cylinder, A, from the area of the horizontal cylinder, which we found 
was .25Tr. The equation for the "back" surface is + « (.25)^ • The vertical 
cylinder is defined by x^ + y^ « (.lO)^, The projection is then half a circle 



of radius .10 (figure 17), and is given by y « /•lO* - x* . 



-.10 .10 



Fig. 17. Projection of the "neck" region 



For a given x, y varies from 0 up to /.lO'^ - x"^ . The limits on x are -.10 to .10. 
The equation for the surface of the back yields the required partial derivatives: 



3z 



f.lO f/. 102-X2 



.1^=0 



-.10 



0 



This strange formula still has visual intutive appeal since the limits on y 
are obtained from the width in the y direction (for a given x) and the limits 
on X are from the length in the x direction. The quantity in brackets accounts 
for the changing height of the surface. The remaining parts of the problem are 
left as an exercise. 
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Error Estimation 



The Integral can be used to develop approximate solutions to certain 
equations or approximate models of given data. The integral provides a qualitative 
error estimate for the approximation. The most well known application is the least 
squares fit of a line through a set of data points. One of the most recent 
applications is the residual norm as an error Indicator for approximate solutions 
to partial differential equations. In each example discussed below, a function is 
integrated over a domain of interest. If this function represents the difference 
between the approximate and true solution, then the value of the integral decreases 
as the approximation improves. The Integral is then minimized to provide the 
"best" approximation. 

Least Squares 

Many experiments produce data as pairs of numbers 



The underlying relationship is often assumed to be linear, that is, the model is 
assumed to be 

y = Ax + B, A,B " constants. 
The points (x^,y^) are usually not collinear (see figure 18) due to experimental 
error, inaccuracies in the model, round-off error during measurement, etc. Thus 
the problem is to choose the constants A,B so the line matches the data points 
as closely as possible. The method of least squares uses the sum of squared 
deviations for the error function, E^ (see figure 18): 



(x^,y^), (X2,y2) 



, . . . , 




n 



(13) 



- I [y(xj - y,]' 



1-1 



where y(x ) « Ax + B. The goal is then to choose A,B so that E^ is minimal. 
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By writing as a function of the parameters A, B, we minimize by setting 
the partial derivatives equal to zero: 

3e2 n - 




y« Ax + B 



deviation* y(xi»)-yi» 



■ Axi» + B - yi 



(xi,yi) 



Fig. 18 Q Least-squares line. 



We then obtain two linear equations in A,B which are easily solved (see problem 5) 
and can be shown to give the line with the least sum of squared deviations (see 
problem 6 ) , 



EKLC 



Norms 

The squared deviation used in the least squares method is an example of a 
norm . A norm of a function f(x), denoted || f(x)|| , is a functional (Pearson 197A, 
p. 9A6) such that, for any scalar a, and any functions f(x), g(x), 

||f(x)|| >0 

||f(x)|| -0 if and only if f(x) = 0 
||f(x)l| - |a| . ||f(x)|| 
llf(x) + g(x)j| < ||f(x)|| + ||g(x)|| . 
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When a curve y « f(x) is approximated by a least-squares straight line 
y * Ax + B over the interval la,b], we choose A, B to minimize the error norm E 
defined by 



[Ax + B - f (x)]^ dx 



(14) 



Note that equation (13) is a discrete analog of equation (14). 

The norm can be used for evaluating the closeness of fit of one curve to 

another, or for obtaining a qualitative estimate of the accuracy of an approximate 

solution to an equation. The L norm of a function f (x) is defined as 

P 



Lp(f(x)) - l|f(x)l| p 



|f(x)|Pdx 



^ a 



i/p 



The least squares norm E is in fact L2(f(x)), since 

1/2 



l|f(x)|| 



if(x)12dx 



and thus 



e2 = [ ||f(x)ll ]2 « [L (f(x))]2 . 

2 2 

Example 8 

An interesting comparison is between the line fitted to the data points of 
example 5 and the line fitted to the smooth exponential curve of eqn. (12), 
rewritten as follows. 



-.05x 



y^ - 7.3(x + 3.5)e 
We can approximate this function with a linear function, 



(15) 



y^ - Ax + B 



3^ 



(16) 
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by choosing A and B to minimize the L2 (least squares) norm of the difference 
between them, * Denote the norm by 



E -Hyi - 7211 - tp(yi-y2) dx]*^ . 



2 '3 

Substituting eqns. (15) and (16) yields 

.30 



E - [| (7.3[xH-3.5]e"-°^'' - Ax - B)2dxl^ . 



(17) 



We now minimize E by taking partial derivatives with respect to A and B (see 
Hertzberg 1977 for a review of differentiation). One first derivative will now be 
calculated, with the remainder of the minimization left as computer exercise 2 . 

First, note that the minimum of E occurs at the same values of A and B which 
minimize E^, since E is positive. The derivative of E^ with respect to A is now 



attempted. First, we write 

9(e") ^_9_ 



3A 



dA 



30 ' 

f (A,B,x)dx 



where f (A,B,x) - (7.3[x+3.5]e~*°^''-Ax - B)^ 



Then the derivative is calculated after the integration is completed, which is not 
a simple task. It may be easier to differentiate under the integral sign first, and 
then integrate the result. 

Theorem (Pearson 197 A, p. 100). 

Let f(x,y) be an integrable function of x for each y, and let 9f/9y be 
continuous over a<.x£b, c <y <_ d. Then 



dy 



f (x,y)dx 



(9f/9y)dx . 
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Thus we evaluate: 



3(E^) 
dA 



30 



3 

rso 



3f A 
3A 



2(7.3[xf3.5]e''°^^ - Ax - B)(-x)dx 



(18) 



which can now be integrated (computer exercise 2 ) • 

No justification has yet been given for the constants in eqn, (15), These, too, 
may be determined by a least squares norm. See computer exercise 1 for details. 

The norm can also be used to estimate qualitatively the relative accuracies of 
approximate solutions to a differential equation. The limits are determined as 
above by the interval in which the accuracy is to be judged. An example is presented 
in the next section. 



DIFFERENTIAL EQUATIONS 

As was previously mentioned, the integral is also used to solve a differential 
equation. In the simplest form, if an equation is written as 



then the solution (general .olution) is 

f 



x(t) 



f (t)dt. 



If we know the value of x at t - 0, the initial condition , then the differential 

equation has a unique solution given by 

ft 



x(t) 



f(T)dT + x(0). 
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Separation of Variables 

When the differential equation also involves x on the right-hand slde^ the 
solution Is not so easily represented • Other techniques must then be employed 
to represent the problem by a set of integrals whose Integrands are functions of 
only the variable of Integration. One such technique Is separation of variables • 
Example 9 

Let us return to example 4. The equation Is 



II « rN(l - N/K). 



(19) 



The simple representation given above falls here: 



N(t) 



rN(l - N/K)dt. 



To Integrate, we need to know N as an explicit function of t, which of course 

we do not know. However, by multiplying both sides of eqn. (19) by the differential 

dt, we obtain 

dN - rN(l - N/K)dt 

and by dividing by N(l - N/K) we successfully separate the variables N, t on either 
side of the equality: 



dN 



N(l-N/K) 



« rdt. 



Integration is now possible since each Integrand Is a function of just the variable 
of Integration, 



N(l-N/K) 



dN 



rdt. 
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How are initial conditions incorporated now? Since the indefinite integral 
produces an integration constant, the constant is determined by the initial 
condition y laaking the general solution unique. 
Example 10 

A recent study on water purification (Dickson, et al, 1977) included the 
mortality of fish due to intermittent chlorination. The following data were 
obtained (Table A) where frequency is the number of doses per 24 hours. 

Table A Fish Mortality Due to Intermittent Chlorination 



Mortality 


Frequency 


Duration (hr) 


.01 


1 


1.5 




2 


0.7 




4 


0.3 


.10 


1 


3.2 




2 


1.4 




4 


0.7 


.20 


1 


4.3 




2 


2.0 




4 


0.9 


.50 


1 


7.5 




2 


3.4 




4 


1.6 



A mathematical model has been devised which shows the effect of duration (T) 
and frequency (F) on mortality (M) using two differential equations: 

3M/3T a^T^'^ , (2 0) 

aM/3F - a2F^*^ (21) 

^1 
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where a^^ is a function of F and is a function of T. Another model was derived 
by curve fitting, using Table A and is quite different. For a given M, 

And) - b^ - b2 An F (22) 

where it is now assumed that b^ is a fxinction of M: 

b^ - 2.16 + 0.40 £n M. 

Can both of these models hold? Are they consistent? We can answer this by 
integrating the first model to derive a single expression for M as a function of 
F, T. Each equation (20, 21) involves a single partial derivative, so that one of 
the variables can be assumed constant. By holding F constant, we integrate eqn. (20) 
to obtain 

M - (a^/l.5)T^-^ + 

which is more precisely written, with a^ as a function of F, as 

a (F) 

M = ^75-r-^ + C^ . (23) 
Similarly, we integrate eqn. (21) with respect to F, holding T constant. 



M 



a2F^-^dF - ia^/l.S)!^-^ + 



a2<T) 2 8 



Thus 



1.8 2 



a3^(F) 2 8 
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We asstime no observable mortality if no chlorlnatlon occurs, so that. If 
F - 0 or T - 0, then M - 0. Thus - C2 - 0. 

,2.5 ^2<^> ^2.8 

Therefore, we must conclude that 

a^(F) - 1.5C I^*^ 
a2(T) - 1.8C T^*^ 

vbeire C is some proportionality constant. From eqn. (23) we have 

M - C F^- V*^ (24) 

It is now straightforward to manipulate eqn. (22) to look like eqn, (24). 

Jin T - (2.16 + 0.40 iln M) - F 
Jin T + b^ iln F - 2.16 



0.40 

b. 



- iln M 



iln T + iln(F h - ilnCe^'^^) _ . 
0.40 

^2 

2.5 iln(TF ^/8.7) - iln M 

Jln[(TF^^/8.7)^*^] - iln M 
2 S 2.5b2 2 5 



.2.5 



Thus, bj - 2.8/2.5 - 1.1, C - 1/(8.7) * - .0045 and the models agree. 
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Residual Norm 



Often a differential equation cannot be solved by any standard methods and 



an approximate solution is developed. In some cases, bounds on the error in the 



approximation can be derived. Too often, however, calculating the error bounds 
is as complicated as determining the solution itself. A recent idea is to use 
the least-squares norm (or occasionally some other norm) as a gauge of the ac- 
curacy of the approximation. With only the differential equation at hand, a 
new type of norm, the residual norm, is used. 

Consider a differential equation written as 



with initial condition x(0) « x^. If the equation is rewritten as 

dx/dt - f(x,t) « 0 

then we obtain an expression which vanishes when the exact solution is used. 
In operator notation the equation is 

L(x,t) » 0 

where L is the differential operator defined by 

L(x,t) « dx/dt - f (x,t) 
If L operates on an approximate solution, x(t), then L(x) ^ 0, and its value is 
the residual . Yet we would expect L(x) to be small if x is a good approximation. 
The final step is to decide on some interval over which the approximation's accuracy 
is to be judged, say ta,b]. Then the least-squares residual norm for the approxi- 
mate solution x(t) is 





a 
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Note that || R || ^ - 0 if the exact solution is used, since L is then 0. This 
limiting behavior suggests the residual norm as a qualitative indicator of the 
accuracy of an approximate solution. No theory exists, however, for estimating 
the actual error from the value of the residual norm. 
Example 11 

A deer population is threatened by severe storms as well as reduced grazing 
area due to a small fire. Consequently, a moratorium on deer hunting is imposed 
until the deer population returns to 60% of the carrying capacity for the area. 
How long should the moratorium be imposed? 

Assumptions are now presented which allow us to model the population dynamics 
of the deer and estimate the duration of the moratoriimi. Let the logistic growth 
model be used: 

^ - Rn(l-n/K) . 
at 

Assume that the current population size is 400, the carrying capacity (K) is 
2000, the intrinsic growth rate (R) is 1.0 and that the approximate solutions 
should be fairly accurate for the first three years. Assume also that the 
moratorium is not obeyed immediately, so that hunting pressure gradually tapers 
off. Let the growth model with the "harvesting" (i.e. hunting) term be 

4^ o Rn(l-n/K) - e"^n , (25) 
dt 

with initial condition n(o)-n . 

o 

Assume also that you do not have access to either a programmable calculator 
or a computer. This equation cannot be solved by separation of variables (try 
it), so some approximation must be made. R:om previous studies, you learn that 
the following two models have been used fairly successfully to model short term 
population growth. 
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n^(t) - (K-n^) d-e"^^/^) + (26) 

HjCt) - n^e (27) 

Note that n^(t) tapers off at the carrying capacity, but n2(t) is exponen- 
tial growth for n^aO; both are properties the exact solution must possess. Which 
approximation is better? We compare n^ and n2 by evaluating their least-squares 
residual 'norms. Rewrite eqn. (25) as L(n) ■ 0, where L(n) represents the differ- 
ential operator acting on n: 

L(n) « dn/dt - Rn(l-n/K) + e"^n - 0 
We expect a good approximate solution, n , to give 
L(n^) « 0 , 

and thus its residual norm should also be small: 

T 



lUCn^)!! 2 =[| a(n^))2 dt]^ * 0 . 



'o 

Intuitively, the smaller || R || ^ is, the better the approximation should be. For 
this example, assume the following values: 

K = 2000 R - 1 

n » 400 T = 3 . 

o 

The integrands in ||R(n^)||^ and ||R(n2)||2 are now calculated, from eqn. 
(26) we obtain, after a few steps, 

[L(n^)]2 = (2000 e"*^- 1600 e"^'^*^ + 1280 e"*^^ 



-2880 e"*^*^ + 1600)2 ^ ^^S) 



Irom eqn. (27) we obtain 



[L(n„)]2 - [ (400/3) e*=/^ - (1-400 e^^^J 



'2' J - I v""'-/ - - vx-^ww c /2000 

-e"*=) 400 e*'^^]^ . (29) 
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After multiplying out the square in each equation and collecting terms, we can 
write both (28) and (29) as sums of exponential functions, and thus both can be 
integrated by hand easily (but tediously). The results are, upon integrating 
over the interval [0,3]: 

II R(nj^)||^ - 661.8 

II Un^)]]^ « 172.3 . 

Figure (19) shows the accuracy of the two approximate solutions as compared 
to a numerical solution to the original differential equation (25). It appears 
that, for the chosen values of R,K and n^, the approximation n2(t) is better, 
and the residual norm for n2 supports this evaluation. The 60% of carrying 
capacity, i.e. 1200, is reached in 3.30 years using n2(t), and in 3.08 years 
according to the numerical solution n(t). It is interesting that, without the 
term for hunting (curve n^(t) in fig. 19), the population reaches the 60% level 
in only 1.79 years. 
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PROBLEM SET 

G eneral 

1. An application of both the definite integral an^ exponential growth is a 

model for the demand for natural resources. Assume P(t) represents the rate 

of use of a resource at time t ^ 0, and represents the rate of consumption 

at time t « 0. The exponential model then gives 

P(t) « P e^^ . 
o 

The constant k can be defined as the rate of increase in the use of this re- 
source. Let A(T) represent the amount of resource used during the interval 
[0,T], T > 0. 

a. Write the definite integral representing A(T) i-n terms of a general 
function P(t). 

b. Evaluate the integral for general T >^ 0 by substituting the above 
exponential function for P(t). 

c. In 1973 (say t « 0), the world use of copper was est ted to be 
6.99 X 10^ kg, and the demand was increasing at an exponential rate 
of 8% per year. How many tons of copper will theui be used from 1973 
to 1983? 

(Hint: First find P , k, T and then use your answer to (b).) 

o 

d. If the rate of growth in demand for coppc"*^ remains? at :jX and no new 
reserves are discovered, when will the world su; , ^ - copper be 
exhausted? Assinne the reserves in 1973 are 336 x 10^ kg, and no 
recycling occurs. 

is 
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2. live thousand trout, each weighing 80 gms, are planted in a lake. The 
population size, N(t), declines exponentially according to the equation 

N(t) - SOOOe^'^^ 

where t is measured in months. Each fish grows according to the formula 

W(t) « 10000 (1 - .8e""°^^)^ 
where W is the weight in grams. 

a. Check that both N(t) and W(t) satisfy the given initial conditions. 

b. Write the formula for the biomass (total weight of all the planted 
trout) in the lake as a function of t, and calculate the biomass 
at 12 months. 

c. What is the average biomass for the first 12 months? What was the 
initial biomass? Sketch hew you think the graph of biomass versus 
time would appear. Mark on the vertical axis the initial, final 
and average biomass values. 

3. Assume that a radioactive element disintegrates so that the number (N) of 
atoms present at a given time (t) decreases at a rate proportional to N itself. 
Let k be the constant of proportionality and be the number of atoms present 
at t «= 0 (thus N(0) = N^). 

a. With the derivative representing the rate of change of N, use the 
above defined parameters and variables to write an equation for the 
rate of decrease in the number of radioactive atoms with passing time. 

b. Solve this different i.a.l equation for N(t) so that the solution satis- 
fies the initial condil';;.Dn, N(0) = N^. (Hint: first separate var- 
iables. Do not forse^: the integration constant I) 
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c. Determine the "half -life" : the time (t^^) such that 

N(t-) » .5 N 
1 o 

At t » t^f half of the original number of radioactive atoms have 
disintegrated (decayed) • 

d. We now make certain assumptions which will allow us to estimate the 
age of a piece of wood by radiocarbon dating: 

i) living organic material contains the carbon isotopes C^^ and 
C^^ in a fixed proportion independent of time, 
ii) the C^^ atoms do not disintegrate, 
iii) the C^** atoms disintegrate radioactively with a half-life of 
5568 years, and 

iv) the C^^ and C^** atoms are not replaced once the organism dies, 
even though the C^** atoms are being lost through disintegration. 
A sample of wood in an American Indian cliff dwelling is measured to 
have 37% of the C^** isotope expected in living wood. Estimate the age 
of the sample. 



50 



46 



4. The survival and health of inter tidal organisms often depend on the amount 
of time they are exposed, i.e. not covered by sea water. A given position on 
the beach can be identified by its tidal height:: the height of the tide when 
the water's edge Just touches that spot. A model of tidal height as a function 
of time can then be used to approximate the length of time between successive 
high and low tides (or vice versa) during which a position on the beach is ex- 
posed. 

The following data were taken from tide tables for Port Townsend, Washington 
for the day of maximum tide difference for daylight low tides. 
Time Height (m) * 

4:35 2.47 
11:49 -0.82 
19:40 2.77 
For each of the two times intervals 4:35-11:49, 11:49-19:40, the model uses the 

following cosine function fitted to the tide data: 

f^fm NN ^ « ti high tide 4- low tide 
H « a cos (b(T-Tj^)) + H^, » — ^ 2 • 

a. Sketch a graph of the cosine function with the maximum at the high 

tide level and a minimum at the low tide level for the first time 

interval. What do the constants a, b, T^^ represent? Calculate these 

constants and H for each of the two time intervals, 
a 

b. The constant H seems to be the average tide height. Verify this for 

a 

the first time interval by evaluating the aJ>propriate definite integral. 
Prove that the average height H for the two intervals combined can be 
written as the weighted average of the average heights H^^, H2 for both 
sub intervals. That is, show that (using decimal hours) 
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H 



15.09 



19.67 



■ Bit -1 iHI »2 



4.58 



c. Ibr a point on the beach at the 0.0 meter tide level, how long is 
it exposed between these successive high tides? (Hint: find the 
inverse function T « g(H), but be careful with the second interval.) 



Least Squares 

5. Write equation (13) using y(x^) « Ax^ + B. Evaluate the two first partial 
derivatives and, by setting both equal to zero, obtain two linear equations in 
A and B. 

6. a. Ibr the following data, derive the "least squares" line and prove 

that for your choices for A and B, the sum of squares is at a true 
minimum: 

f(x,y) is at a minimum at x^, y^ if, for x«x^, y^y^ » 
i) 9f/9x » 9f/9y « 0 , 
ii) 02f/9x2)(92f/9y2) - (92f/9x9y)2 > 0 , and 
iii) 92f/9x2 > 0 . 
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Data froQ flg» 8 In text , 

X y 

3.0 39.0 
5.3 55.0 
9.6 57.5 
18.0 60.0 
30.0 57.0 
Now determine a **best'^ exponential fit to the above data. FLrst 

transform the desired function 

A Sx 
y « Axe 

Into a linear function (with respect to x) 

h(x,y) » f (A,B) + g(A,B)x. 
A neater format Is 

h « f + gx. 

Make a table of values for x and h. Determine the least squares 
line through this new data set and. In the process, calculate f , 
and thus A, B. 
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Answers to the Problem Set 



a. 



A(T) represents the total change In the resource < 
rate of change, P(t), is Integrated: 

T 

A(T) - P(t)dt. 



Thus the 



b. With P(t) - Pj^e*'^, 



T 

f 



A(T) 



o 



kt 



P^(l/k)e 



(P^/k) (e*"^-!) 



c. The initial use rate is P , the rate of increase is k, and the time 

o * 

interval ends at 1983-1973 = 10 years. Thus 

9 

P = 6.99 X 10 kg 
k » .08 

T = 10.0 years, 

and from (b) , 

A(10) - (6.99 X l0V.08)(e^*"^^^^°^-l) 
- 107. X 10' kg 
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We solve for T. The total change equals the current 
(at t ■ 0); thus the parameters are 
Pq - 6.99 X 10^ kg 
k - .08 

A(T) « 336 X 10^ kg 
and from (b), 



336. X 10' 



(6.99 X 10«)e-°^*'dt 



= (6.99 X l0V.08)(e-°^^-l) 
3.85 = e-°8T . 1 
4.85 = e-OS^ 
in 4.85 - .08 T 

T = (in 4. 84)/. 08 >= 19.7 years. 

The initial conditions are N(0) - 5000, W(0) - 80. 
N(0) = 5000e~'^^°^ - 5000(1) » 5000 
W(0) - 10000(1 - .8e~-°^^°^)' 
« 10000(1 - .8)^ 
= 10000 (.2)^ « 80. 
Let B(t) represent the total biomass. 

B(t) « N(t)«W(t) 
With t « 12, 

B(12) = 5000e--^^^2\l0000)(l.-.8e--°^^^2^)3 
= 21.88 X 10^ g. 
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Average blomass B(12) uses the definite integral 

12 



¥(12) = (1/12) 



B(t)dt. 



Substitution for B(t) gives 

12 

B(12) = (1/12) I 5.0 X 10'e"-^^(l. - .8e"-°^^)3dt 

0 

12 



4.17 X 10' 



-.5t,, o -.05t.3 , 
e (1. - .8e )' dt. 



Integrate by parts, where 



o -.05t.3 

u « (1 - .8e ) 
dv - e dt. 



12 



e--5t(l - .8e--°^^)^dt 



^ (1 - .8e-°^^)3 



-.5 

12 



12 



-.5t 



(3)(1 - .8e"-°^^)2(.04)e"-°-^dt 



12 



- (-.00496) (.177) + .016 + .24 



e (l-.8e 
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2. 



c. (contd.) 



Integrate by parts again, with 

/n Q -.05t.2 

u = (1 - .8e 
dv » e '"^^dt 

An intermediate stage is (two steps skipped) > with the left side included, 



12 

r 



£ (1 - .be ) dt 



.0324 + .0349 



0 
2 



.0324 + .0349 



-.60t o -'^St,^ 
e - .8e dt 



= .0324 + .0349(.4353) 
= .0476 

B(12) - (4.17 X 10^) (.0476) - 1.98 x lO^g . 
You might assume the function is monotone decreasing. The rapid 
initial weight t^ain, however, causes a slight initial rise, as 
shown below. The initial biomass is the product of initial weight 
and initial population: 

B(0) « W(0)N(0) = (5000) (80) = 400xlO' gm = 400 kg. 
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600r 




4 6 8 10 12 
TIME (MONTHS) 
FIO. Ifl. TROUT BIOMflSS MODEL 



The rate of change of K is a decrease (I.e. negative) and proportional 
to N. For k > 0, we have 
dN 



dt 



-kN 



Separate variables. 



dN 



dN 
N 



-kdt 



An N - -kt + C, C - integration constant 

-kt C 
N " e e 



But N(0) -= Nq. Thus 
N(0) « Nq = e°e^ 

and the solution Is 

N(t) - N^e"*'*^ . 
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3. c. Substitute Into the solution equation. 

N(tj_) - N^e 

-kt 

Jin (.5) - -kt^ 

- -an(.5))/k 
-= Un 2)/k 

d. Let N(t) be the number of C^** atoms In the wood. Then 
N(t)/NQ -= .87 . 
The half-life of 5568 gives 
Un 2)/k «= 5568 

k «= (£n 2)/5568 - .00012A5. 
The equation for N(t) then gives 
N(t)/NQ - e"^^ 

.87 - exp(-.00012A5t) 
iln(,87) - -.00012A5t 
t ■ 1119 years 
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■1' ' ' ' 1 1 1 I 

4 6 8 10 12 14 16 18 20 

TIME (HOURS) 

FIO. 2fl. COSINE FITTED TO TIDE DflTfl. 



a » half the total change in height (amplitude of cosine function) 

- (high - low)/2 

b = 27r/p where p « period « 2|t^^ , - T, | 

' high low' 

= "phase" lag. The simple cosine function begins at t"0, 

so the lag is the time of high tide. 

First interval: 

a « 1.645m., b - 27r/(2|4.58-11.82|) - 0.434 (decimal hours) 

T, - 4.58, H « 0.825m. 
h a 

Second interval: 

a - 1.795m., b - 27r/(2 |l9.67-ll. 82 1 ) - 0.400 

T, - 19.67, H - 0.975m. 
n a 
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b. H average height, first Interval "7.24 hours 

11.82 



7.24 



a cos (bt - bT. ) + H dt 
n a 



4.58 



(-a/b)(sin(11.82b - 4.58b) - sin (4. 58b - 4.58b)) 



7.24 ^ a 

a 



7.24b 
+ H 



(H ) (11.82 - 4.58) 
(sin(7.24(2TT/2(7.24))) - 0) 



a 



" H , since sin it « 0. 
a 



For tlie two intervals combined, we must weight the average height 
in each interval by the length of the time Interval. 

19.67 11.82 19.67 

H(t)dt 

4.58 



H = 



15.09 



15.09 



H(t)dt + 



15.09 



H(t)dt 



4.58 



11.82 



11.82 



19.67 



7.24 1 



15.09 7.24 



H(t)dt+^f ^ 



H(t)dt 



4.58 



11.82 



7.24 +745- 



15.09 1 15.09 2 



I] 

6 
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c. We need the times of day when the 0.0 tide level is reached. 




"4 6 8 10 12 14 16 18 20 

TIME (HOURS) 

. k EXPOSED*! 
FIG. 3fl. EXPOSURE TIME AT O.OM TIDE HEIGHT. 

The easiest way is to invert the function H(t) for each subinterval. 

H - a cos(b(T - T, )) + H 
n a 

(H - H )/a - cos(b(T - T, )) 
a n 

cos"-^((H - H )/a) « b(T - T, ) 
a n 

(l/b)cos"-^((H - H )/a) + T, - T 
a n 

Now substitute the appropriate constants, and 0 for H. 
ilrst interval: 
T 



cos"-"- (-.825/1. 6m5) + '4.58 



.A3A 

9.A1 « 9:2A a.m 



Second interval: 
T = 



•'■ cos"^(-. 975/1. 795) + 19.67 



.AGO 
25.033. 



This answer is obviously v;ron6 since the 0.0 height must be attained 
some time between low tide (11.82 hours) and the next high tide 
(19.67 hours). The fallacy is in treating the "arccos x" as an 
inverse function. It is one-to-one only when the domain of "cos x" 
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c. (contd.) 

is restricted to [0,1:]; I.e., arcos x always has [0,7r] as Its range. 
Thus the first term 

cos"^(-. 975/1. 795) 

correctly tells how far from the maximum (T^^) is the desired time T, 

but not whether T lies above or below T. . With the second subinterval, 

n 

T is obviously below T. and hence we subtract from T. : 

n n 

T = 19.67 - -—^ cos"-'-(-. 975/1. 795) 
= 14.31 hrs. - 14: 18 . 



n 

E = ^ [Ax,+B - Y ]2 
i=l ^ 
n 



3E 
3A 



i«l 



Set equal to 0: 



n 

I 2[Ax +B - y ]x - 0 



I [Ax.+B - y ]x. - 0 (a) 
1=1 ^ 



Equate to 0 and divide by 2: 



n 



I [Ax +B - y.] = 0 (b) 
1=1 
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5. (contd.) 

Thus we have two equations which are linear in A and B. This 
is more easily seen if we rework them. Rrom (a) 



n n n 

.2y 



A( Z + B( 2 X.) - ( 2 x.y.) - 0 . (a*) 
i"l ^ i-1 ^ i-1 ^ ^ 

And (b) becomes 

n n 

A( 2 X ) + n B - ( Z y ) - 0 . (b*) 
i=l i-1 ^ 

a. Eq. (a*) is: 

A(1353.25) + B(65.90) - (3750.50) - 0. 

Eq. (b*) is: 

A(65.90) + 5B - (268.50) = 0. 

We solve for A and B. 

B = (268.50 - 65.90A)/5. 

A = (3750.50 - 65.90B)/(1353.25) 

" 2.771 - 2.615 + 0.642A - 0.156 + 0.642A 

Thus 

A ■= 0.156/(1. - .642) = .436 
B " (268.50 - 65.90(.436))/5. = 47.95. 
To prove that these values for A,B do give a minimum sum of 
squares, we must show that, for A = .44, B ■= 47.95, 
i) 9f/9A - 9f/9B = 0 
ii) (92f/9A2)(92f/9B2) - (92f/9A9B)2 > o 
iii) 9^f/9A^ > 0, 
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a. (contd.) 



where f is the sum of squares 



f - Z(y - Ax - B)2 . 

i 



We have 



3f/9A - -2 Z(y - Ax - B)x 

i 

3f/3B - -2 Z(y - Ax - B) . 

i 

Since A,B were calculated by setting these derivatives to zero, we 
do not need to check (i) . 

S^f/SA^ = 2 Z x^, 32^/332 , 2n 
i 

3^f/3A9B = 2 Z X . 

i 

Substituting into (ii) gives 

(2 Zx^)(2n) - 4(S x,)^ - 20 SxJ - 4(Zx 
i^ i i 

= 20(1353.25) - 4(65.90)^ 

« 9693.76 > 0 . 

Note that (iii) is obviously satisfied so we have a miniraum for f. 



b. Begin with 

A Bx 
y = Axe 

Divide by x to isolate the exponential function. 
Bx 



y/x «! Ae 
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b. (contd.) 

Now take logarithms of both sides. 

AnCy/x) - £n A + Bx 
Thus, In the form h » f + gx, we have 

h - iln(y/x), 

f « £n A, 

g - B. 

The transformed table is below. 

X h 
3.0 2.565 
5.3 2.340 
9.6 1.790 
18.0 1.204 
30.0 0.642 
The least squares calculations give the transformed equation 
h » 2.642 - 0.071X 

and thus 

- , n/i -.071x 
y « 14.041X e 



c. The following graph (fig. 4A) compares the least squares line 
and exponential curve of parts (a) and (b) with the data 
and the exponential curve in the text (eqn. 12). The fit 
of eqn. (12) seems superior to that of the curve from (b). 
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(co: tel.) 

However, the fit of 

h - 2.642 - 0.071X 
to the transformed data is quite good, as shown by fig. 5A. 
Thus we conclude that a good least squares fit to transformed 
data does not necessarily imply a good curve fit to the 



original data 
75r 



least squares exponential fit 
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5 10 15 20 
OXYGEN CONSUMPTION 
FI0.4fl. COMPARISON OF FITTED GROWTH FUNCTIONS. 




n' ■ I I I I I 

0 5 10 15 20 25 30 
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COMPUTER EXERCISES 



Program INTG Is both a calculation and graphics routine « Simple Integrals 
are not calculated. The solutions are predetermined and are plotted or tabulated 
as a check on the user's own calculations. Integrals which have no known closed 
form solution are calculated numerically, using an adaptation of Newton-Raphson 
or Romberg schemes. 

The specific user options and restrictions are detailed In the section of 
this monograph titled "User's Guide for Program INTG." This program Is designed for 
the line printer graphics subroutine PRNT3D, whose operation Is described In Gales 
(1979). In each exercise, the user must Input certain paramete,r values. Some 
parameters control the behavior of the functions. By repeating an exercise with 
different parameter values, the user can gain a better Intuitive understanding of 
the degree of dependence a given function or method has on the various parameters. 

1. The function used In eqns. (12) and (15) was determined by trlal-and -error 
using a calculator. The values of the three function parameters were changed until 
a combination gave an acceptable visual fit to the data points. A more consistent 
approach Is to use least squares. The function Is 

y - C(x+A)e"^^ 

so the goal Is to choose the parameters A, B, C to minimize the error function 
E, given by 

E - Z (y . - C(x +A)e"^''l)^ 
1«1 ^ 

The difficulty with the use of a nonlinear least squares approach is that 
the equations for the unknown parameters often must be solved numerically. 
Program INTG uses a Newton-Raphson iterative scheme to evalr'^te the parameters. 
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Starting values are inputted to the routine which iterates until the values are 
relatively constant between successive iterations. The problem is that the 
Iterates may never converge, or may converge to a local minimum for E which is 
much greater than the absolute minimum, and which determines values for A, B and 
C that are far r*>moved from thtlr "optimal" values • 

When the approximation is used only for a curve f say for estimating 
derivatives, the disparity between the compuced vs. optimal values for the 
parameters is of little consequence. However-, when the function is used to 
suggest a biological mechanism, i.e. each paramete*' has a biological interpre- 
tation, the differences can lead to erroneous explanations of the data. 

program INTG illustrates the sensitivity of the nonlinear least squares 
method to the initial parameter estimates. Use OPTION-1 and try the parameter 
values used in eqn. (12) as your initial estimates. Your output will include 
a table of values for the iterations on A, B and C, to demonstrate the degree 
of convergence. The output will also include the value of the residual sum of 
squares and a graph of your approximate function and the data points. Repeat 
the exercise with different starting values and try no achieve a lower sum of 
squares. Be sure to try (A, B, C) ^ (5.0, .05, 5.0) and (6.0, .05, 5.0). To 
what value do the iterations converge? What 1p the graph of the second set? 
My minimum sum of squares was 52.0. Can you find s parameter set (A, B, C) which 
gives a lower sum of squares? 

2. Another way to use least-squares is to approximate one function by another. 

Here we approximate the curve 

c . r n\ '0.045X (30) 
y « 6.1 (x+5.0)e 

by a straight line 

y « Ax + B • 
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The goal Is not to minimize a sum of squares, but to minimize an integral, also 
denoted E, called the "least-squares norm" (see example 8, p. 33). Whereas 
example 8 used the function 

y - 7.3 (x+3.5)e"'°^^ , 
we will use the above curve, since it gives a better fit to the experimental 
data (see exercise 1, above). 

a. First, finish example 8 using the above function (eqn. 30), i.e. evaluate 
both 9(E^)/9A and 9(E^)/9B, then set both equal to zero to obtain two equations 
in A and B. The goal is to minimize: 

.30 



E^ 



[6.1 (x+5.0)e^-°^^^-Ax-B]2 dx (31) 
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3 

(Hint: when evaluating 9E^/9A and 9E^/9B, differentiate under the integral sign.) 

b. Now, before solving these two equations for A and B, look at the graph 
from exercise 2 or even figure 7 and choose values for the constants A and B 
which you believe will describe a line which is a fair approximation to the 
curve. Now use program INTG with 0PTI0N=2. Your input includes values for A 
and B. This option causes the least-squares integral to be calculated using a 
Romberg numerical Integration scheme. The output is a table of the iterations 
of the least-squares integral (to illustrate the convergence of the numerical 
routine) and a graph of the data points, the function of eqn. (30) and your linear 
approximation. You may wish to repeat the 3x&rcise with several choices for A 
and E in an attempt to improve the fit. Now solve your equations for A and B 
and input them into program INTG. How close were your previous guesses? Does 
a dramatic reduction in the integral (31) correspond to a much Improved fit? 

3. The stomates in leaves are small pores which permit the exchange of gases 
and water vapor. The stomates change shape from a near circle to a slit In order 
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to adjust this exchange in response to varying environmental conditions. Each 
stomate is roughly elliptical in shape with constant perimeter throughout the 
shape change. The area of the opening is then given by 



A 2b 
A » — 

a 



a 

(a^-x^)^ dx 

-a 



where 2a, 2b are the lengths of the major and minor axes, respectively, of the 
ellipse. 

a. First evaluate this integral using either of two trigonometric substi- 
^tutions: let x » a cosG or let x « a sinG. Be 6i.re to check the limits 
on e and the sign of each trig function in the transformed integrand. Your re- 
sult must be non-negative since the area is non-negative. For a fixed perimeter, 
say 35 y, the area should be maximal when the hole is a circle, i.e. a»b. Demon- 
strate this by evaluating the area when a«b«5.57, and when a-5, b»6.09. (The 
perimeter formula is P «= 27r/a^+b^)/2 .) 

b. Use your area formula (the answer to part a) to investigate the dilution 
of the sun's energy flux due to the angle at which the rays strike the earth's 
surface. The dilution is caused by a fixed amount of solar energy being spread 
over a larger area (as the rays become more horizontal) and thus the energy per 
unit area decreases as compared to the energy density of vertical rays. 

The situation is pictured below, giving a top and side view for two in- 
clination angles. Th6 sunlight passes through a circular hole in an opaque sheet 
and strikes the earth at an angle R from the perpendicular. 
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sunlight 




^sunlight 



side 
view 



'ground surface 



top 
view 





The Illuminated area changes from a circle (R"0) to an elongated ellipse 
(as R increases). For a given diameter B, write the illuminated area A as a func- 
tion of the angle of inclination R« Check your result by choosing values for 
B and R and calculating the dilution D, which equals the circle's area divided 
by the ellipse's area. Now use program INTG to check your calculations, using 
0PTI0N«3. The input includes B and R. Your constraints are 0 ^ B ^ 5, 
0 ^ R ^ 1.57 radians. The output includes the dilution coefficient and a graph 
of the elliptical region. Does D « .6 mean more dilution (less energy density) 
than D « .3 or vice versa? The inclination angle in Minnesota changes from about 
20** in summer to 65^ in winter. What is the-resulting change in energy density, 
measured by D? 
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SOLUTIONS TO THE COMPUTER EXERCISES 

1. The starting parameter set (A, B, C) • (5.0, .05, 5.0) caxxses 
convergence to the values (5.0, .045, 6.2), with sum of squares equal to 
52.2. The starting set (A, B, C) - (6.0, .05, 5.0) iterates (50 times) 
to the values (12535, ,005, .004), with sum of squares equal to 454. The 
graph for the latter set is essentially a straight line, and a terribly 
poor fit. 

2. a. First partial derivatives: 

3 (6.1(x 5.0)e-°^5- _ ^ _ gja^^ 
dA 3^ OA 

» 2 f (6.1(x + 5.0)e '^^^^ - Ax - B)(-x)dx 

3 

Equating to zero and dividing by 2 gives 

.3 0 ^ _ n45v o 

0 - - f 6.1(x2 + 5x)e -^^^^dx + / Ax^ + Bx dx . 

3 

Successive integration by parts (for the first integral) and direct 
integration (for the second integral) yield the first equation: 
0 - 8991. OA + 445. 5B - 25984.9 

. r^it.lOi -h 5.0)e-°^5^. Ax - B)^dx 

o t^^ ic w ^ ^ r^^ -•045x- Ax - B)dx 
» -2 J (6.1(x + 5.0)e 

Equating to zero, dividing by 2 and solving yields the second equation: 
0 « 445. 5A + 27. OB - 1568.8 . 
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b. The solutions are 

A - .061, B - 57.1 . 
The least squares integral using the approximation Y « .061x +57.1 
is 5A6.452. 

3. a. Let x - a sin 0. The limits of integration are then changed as 
follows: 

~a 4> X 4 B, 

-a ^ a sin 0 < a 

-1 < sin 0 « 1 

There are many options for the range of 0. We first change the integrand 
to see which range is suitable. We know the integrand is non-negative 
since x^ $ a^* Substitution gives 
X « a sin 0 
dx » a cos 0 d 0 

/I^"T1^2jx - /a^ - a^sin 0 a cos 0 d 0 
The lntegr?ad is then non-negative when the positive square root is used 
and wher> cos 0 is non-negative. Thus, we choose the limits of integration 
to be ^2*2 

sin (- I ) « -1, sin ( 5 ) - 1 and 

cos 8 < 0 if - I ^ 0 ^ 5 . 
Ti\t Irwtegral is t-hen 

s; /? ^ TT n _ 

/ /a^- a^sin^0 a cos G d 0 - a^ / /l - sin^0 cos 0 d 0 
-7T/2 ^/'^ 
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Thus 



b. 



Tr/2 



/ cos^© d 0 - f + 

-TT/2 ^ 



G , sin 20 



TT/2 
-TT/2 



1 



- [ TT/4 + (1/4) (0) - (-TT/4) - (1/4) (0)] - TTaV2 
A - ~ fVa^ - dx - Ttab . 



-1 



When a - b - 5.57, A - 97.47 . 

When a » 5., b - 6.09, A - 95.66 . 
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From the right-hand diagram, we have 
cos R « B/x 

and thus 

X » B/cos P . 

The minor axis is then the circle's diameter, B. In the answer to part 
"a", we replace 2a by B/cos R and 2b by B. 

A « Tiab « 7r(B/2cos R)(B/2) 
« TTB^/(4 cos R) . 
Since the area of the circle is 7TB^/4,the dilution coefficient is 



TTB^/CA cos R) 



D « ^^2 /7/ «N - cos R • 
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Thus D - .6 means less dilution than D - .3, and D ■ 1.0 means no 

i 

dilution at all. For Minnesota, the coefficients are 
D - .94 for R - 20° 
D - .42 for R - 65° . 
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USER'S GUIDE FOR PROGRAM IKTG 



Identification 

INTG - A program which determines the values of certain integrals to check 
the user's own calculations or to show the dependence of each solution on 
various parameter values • 

Authors R« Hertzberg, M. Bailey, and L. Anderson, Center for Quantitative 
Science in Forestry, Fisheries and Wildlife, University of Washington, 
Seattle. April 1979. 

Purpose 

Program INTG Is the computer supplement to the module "Calculus- 
Integration" by Richard C. Bertzberg. This module reviews the basic concepts 
of Integral calculus using examples from various areas of ecology, and 
discusses some recent uses of definite integrals in error analysis, namely, 
residual norms. INTG may be applied to three problems: 

1) a numerical, nonlinear least-^squares fit of a function to data points 

using a Newton-<*Raphson iterative scheme, 
ii) a numerical least-squares approximation of an exponential type of 

function by a straight line using a Romberg iterative scheme, and 
ill) the determination of the area of an ellipse, representing the dilution 

of solar energy flux resulting from the rays striking the groxmd surface 

at a nonzero angle of inclination. 

All the problems illustrate the dependence of the solution on 
certain Input parameters. Problems 1 and 2 require initial estimates for 
the function parameters. Problem 3 requires the angle of inclination and 
the cross-sectional diameter of the solar rays. 

Operation 

INTG reads in any number of data sets provided by the liser and 
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■aerates tables and/or printer plots (see Gales, 1978) until it encounters 
mis directivs (see Input ) . 

The operation of INTG is primarily controlled by an input variable 
turned OPTION which may take on three values: OPTION-1 selects a function fit 
to data points, OPTION-2 selects a line fit, and OPTION-3 selects solar 
dilution. For OPTIONS 1 and 2, the' program iterates the numerical approxi- 
mation 50 times, or until the residual becomes too large. For OPTION-3, 
the program evaluates the exact solution. 

Since INTG automatically generates valid default values for each data 
set, the user may override some, all, or none of those listed in the INPUT 
TABLE. Thus, INTG will run correctly without the user specifying any input 
values at all. 
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Program Organization 

The program Is organized according to the following fJ 

' Yes 



lart. 



NODFLT - .T.? 



^ No 



Read default values for 
all variables from the 
built-in default file 



Read in the next 
user-supplied data set 



Output 
appropriate 
error 
messages 



FINIS « -T.? 



Yes 



No 



Check for errors In 
the Input set just read 



Yes 



Were errors found? 



No 



Calculate the appropriate 
Integral; Iterate the 
calculations for 
0PTI0N=1,2. 



Write out X, Y, Z 
coordinates for plots or 
tabular output plus the \ 

titles for the plot. i 



Call the printer plot 
subroutine QQPR3D 
which generates the plot 
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Terminate 
program 
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Input 

-All input is handled by a format-free input package (Gales and 
Anderson, 1978) which permits a user to assign values to variables by a 
^'name»value" convention. No variable need be explicitly assigned by the 
user 9 however, as unassigned variables automatically assume default values • 
The input consists of any number of data sets, each of which is terminated 
by a dollar sign ($). Comments in data sets are enclosed by slashes (/). 
Each data set generates a separate printer plot. 

The input for INTG is divided into three classes: (1) variables 
having mathematical or physical significance: OPTION, A, B, C, DIAM, ANGLE; 
(2) variables which control certain program operations, such as program 
termination or the handling of default input: IPRINT, ECHO, NODFLT, and 
FINIS; and (3) variables which control the printer plots (default values are 
in parentheses" : XMIN (0) , XMAX (0) , YMIN (0) , YMAX (0) , ZMIN (0) , ZMAX (9) , 
XRICH (0), YRICH (0), DFAULT (0), OVPRNT (.F.), AVE (.F.), INT2D (.FO, 
NX (60), NY (45), and ZMAP (0,1,2,3,4,5,6,7,8,9). The variables in the first 
two classes are explained in the following INPUT TABLE whereas the printer 
plot variables are explained in the user's guide for PRNT3D (Gales, 1978). 
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INPUT TABLE 



Name 


Type and 
Dimensions 


Range Limits 


Descriptions 


OPTION 


Integer 


1,2,3 


Number of problem selected: 

1 » function fit to data points 

2 » line fit to a function 

3 * solar dilution 

Default value Is: OPTION « 1 


A,B,C 


Real 


0.1,10.0 


Parameters of fitted function. 
A, B and C are used In 0PTI0N«1; 
A and B only are used In OPTI0N«2; 
and none are used in uriJ.ON»j. 
Default values used are A « 5.0, 
B - 0,04, C - 6.0. 


DIAM 


Real 


1.0,5.0 


Diameter of solar beam used for 

OPTION » 3 only. 

Default value Is: DIAM » 1.0. 


ANGLE 


Real 


0.0,1.570763 


' Angle of Inclination of solar beam 
used for OPTION - 3. 
Default value Is: ANGLE « 0.6. 


IPRINT 


Logical 


• F • f • X • 


A logical value wnicn causes 
the current values for all Input 
variables (default as well as 
current user Input) to be printed. 
Default value Is: IPRINT » .F. 


ECHO 


Logical 


• F« f • T • 


A logical value which causes the 
user's Input to be echoed If 
ECHO ** .T., or suppresses echoing 
If ECHO « .F. 

Default value Is: ECHO » .T. 


NODFLT 


Logical 


• F m f 


1 A logical value wnicn suppresses 
the Input of default values If 
NODFLT » .T. 

Default value Is: NODFLT « .F. 


FINIS 


Logical 


• F • 9 • T • 


A logical value which causes 
program termination If and only 
If FINIS - .T. 

Defaxxlt value Is: FINIS » .F. 
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The last four variables deserve special mention: 

1. The logical variable IPRINT controls the output of all input variables 
which are currently in effect (default values as well tis thor.e specified 
In the current input set). Setting IPRINT-TRUE (or T or displays 
the current values for all input variables; setting IPRINT^FALSE (or F 
or .F.) suppresses the display. 

2. The logical variable ECHO controls the echoing of the input cards. 
Setting ECHO-TRUE causes the subsequent input set to be echoed; 
setting ECHO-FALSE suppresses the echo for the subsequent input set. 

3. The logical variable NODFLT can be used to inhibit the automatic 
assi ^ of default values to input variables. If NODFLT is set 
TRUE 1. >..le current input set, then the current input set is assigned 
default values as usual, but all subsequent input sets merely accumulate 
more input values. In effect, the input values which exist after the 
i-th input set is read, become the default values for the (i+l)-th input 
set. The standard default values may then be restored by setting 
NODFLT-FALSE, but, again, the effects of this change are delved until 
the next input set is read. To a limited extent, NODFLT permits a user 
to set up his own default values and can be very useful for executing 

a number of input sets which differ only in a few parameters. Consider 
the following example in which a user wishes to vary the angle of inclina- 
tion of the solar beam, but keep all other parameters the same: 
/INPUT SET 1: THE FOLLOWING VALUES BECOME THE DE FACTO/ 
/DEFAULTS FOR ALL SUBSEQUENT INPUT SETS; / 

NODFLT - TRUE, OPTION • 3, DIAM - 1.4, ANGLE - 0.2, $ 
/INPUT SET 2: UP THE ANGLE: / 

ANGLE - 0.4, $ 
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/INPUT SET 3: UP IT AGAIN: / 
ANGLE - 0.6, $ 

/INPUT SET 4: YET AGAIN: / 
ANGLE - 0.8, $ 

/INPUT SET 5: NOW STOP / 
FINIS - TRUE, $ 

4. The logical variable FINIS controls program termination. The user must 
add the card: 

FINIS - TRUE, $ 

as the very last input set. If FINIS is not set, the program will termi- 
nate abnormally. 

Output 

For each option, INTO produces one or more tables with convergence 
or supplementary information, and a printer plot (Gales, 1978) which contains 
a title, X and y axis annotation, printer plot lines, and a plot legend 
which helps numerically interpret the plot. The latter is best explained 
by an example. Consider the plot in the first sample run which displays a 
function fit to a set of data points. The function and parameter values are 
specified in the title, as follows: 

NP^TLINEAR LEAST SQUARES FIT OF 

Y = C * (X+A) * EXP (-B*X), (2) 
TO COHO SALMON GROWTH DATA (1) 
A » 5.48 B - 0.04 C ■ 5.67 
The coho growth rates are indicated by the scattered "1" points in the plot, 
and the fitting function by the densely packed "2" points. For example, 
a coho point ("1") near the top middle, of the plot is positioned between 
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1.627 and 2.085 on the x-axis and above 5.962 on the y-axis. Since^the scale 
factor for both axes is E+01 (i.e., lO"*"), the data point represents a growth 
rate of about 60,0 mg/g/day for a dissolved oxygen value of about 18.1 mg/1. 
The last two lines of the legend specify the value of the z-axis levels and 
the number of points mapped to each level. In this case, a large number 
(-9 means >99) were mapped to level 0, five were mapped to level 1 (and show 
up as "l"-s), 64 were mapped to level 2, and so on. 

Restrictions 

The user must restrict all input variables according to the range 
limits listed in the input table, in order to avoid unrealistic physical 
values. 

Error Messages 

There are three types of errors which may occur when attempting to 
execute program INTG: 

1. Syntax errors in the user's input 

2» Range check errors 

3, Plot parameter errors. 
For type 1 and 2 errors, the program flags the error, skips the calculations 
and plotting, and then reads in the next input set. For type 3 errors, the 
program suppres^s plotting, outputs an error message, and reads the next 
data set* For a complete description of type 1 and type 3 error messages 
and actions, refer to the user's guides for the format-free input package 
and printer plot package, respectively. 

Type 2 error messages are of the form: 
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ERROR NO, n IN INTG— 

V " X 

V MUST BE GREATER THAN y 

V MUST BE LESS THAN z 

Where V is the number 2, 3, 4, 5, or 6; v is the name of a variable (A, B, 
C, ANGLE, or DIAM), and x, y, and z are numbers. 

Error message number 1 is slightly different: 

-r-^-r-ERROR NO, 1 IN INTG-r — 
OPTION « X 

OPTION MUST BE 1, 2, OR 3 

Sample Runs 

The annotated listing starting on the next page illustrates the 
control cards and input cards for three sample runs. Each input set must be 
terminated by a $ and generates one or more tables and a printer plot. 
The output appears on the next few pages. 
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XI*NT6#T10#P2. 
ACCOUNT#9XMO8C00# 
COHKENT. 
COMMENT. 
COHHENT. 
COMMENT. 
COMMENT. 
COMMENT. 
COMMENT. 
COMMENT. 
COMMENT. 

ATTACH^ BINTG>I0-BINTG. 
ATTACH#BFF# ID-BFF. 
ATTACH* B PR I0-BPR30. 
COMMENT. 



INPUT SETUP FOR PRJGRAM -INTG- 



********* 

* THE FIRST CARD SPECIFIES THE JDB NAME* 

♦ * -XINT6-* THE CENTRAL PROCESSOR TIME* 

♦ , 10 SECONDS* AND THE PRIORITY (2) * 

* THE MEMORY IS 6DDDD OCTAL bY DEFAULT.* 



* 
* 
* 
* 
* 
* 



COMMENT 
CU.iMENT 
COMMENT 
COMMENT 
COMMENT 
COMMENT 
COMMENT 
COMMENT. *♦*♦ 
COMMENT. 
LOAD^BlNTG. 
LOAO#tlFf . 
LGA0#BPi^3D. 
EXECUTE* INTG. 
COMMENT. 
COMMENT. *♦*♦ 
COMMENT. 
COMMENT. 
COMMENT . 

COMMENT. 

COMMENT. 

♦EQR 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 



THE ABOVE CARDS ATTACH PROoRAM -BINTG* 
(THE BINARY OF PROGRAM -INT6-)# AND * 
THE SUPPOJkT ROUTINES -BFf- (FitEE * 
FURM INPUT)* AND^-BPR3D- (THE PRiNTtK* 
PLOT PACKAGE PRNT3D). THEY ARE ALL * 
IN BINARY FORM ' * 



♦ THE ABOVE CARDS LOAD ThE 

♦ THE SUPPORT ROUTINES AND 

♦ TO -INTG- FOR EXECUTION 



PROGRAM AND* 
PASS CONTRGc 



THE FOLLOWING DTFaULT VALUES ARE ASSUMED UNLESS 
OVERRIDEN tiY EXPLICIT INPUT VALUtSl 



•HkUGRAM PARAME TERS* 
UPTION « l9 

A • 5.0* B • 0.0<i# C • 6.0# 

ANGLE • 0.6# -DIAM ■ 1.0* 
ECHO - TRUE* IPRINT -F AL SE # NODFL T 
FINIS • FALSE* 
•♦***PLOT PARAMETERS***** 

NX • 60* NY • <»5* ZMAP 

XMIN • 0* XMAX • 0* YMIN 

iMlN • 0* ZMAX - 9/ XRICH 

DEFAULT • 0/ OVPRNT • F* AVE • 



■ TRUE* 



» 0* 1*2*3* <t* 5*6* i')' :i>V» 



0* 

a 0* 
F* 



YMAX ■ 

YRICri 

iNT2i) 



0* 
I 0* 
• F, 



/*«**»«*«*****«***«******RUN 1< 

/ 

A • 5.<i8* b • O.Ci* 

/ 
/ 

/***««««««*«**«**********RUN 2' 



5.fa7* 



BR. 
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/ f 
OPTION • 2> A ■ O.TO# 51. Of S 

; 

/ / 

OPTION • It ANGLE • 1 .0# DIAM • l.Oi 

XRICri • .05#YRICH • .05» S 
/- ' 
I ' 
/ 4****^****** Q? PR OGRAH****** 

/ ' 
FINIS • TRUEfi 

♦ EOR 



( 
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QGRAH -INTG- REAOf FOR INPUT 



. THE FULlOWlNG DEFAULT VALUES ARE ASSUMED UNLESS 
OVERRIDEN 8Y EXPLICIT INPUT VaLUESi 

♦♦♦♦♦PROGRAM PARAMETERS'***^ 
OPTION • I, 

A ■ 5,0# B B 0.0<»* C ■ 6.0# 

ANGLE ■ 0.6» DIAM » 1,0/ 

ECHO « TRUE* IPRINT -FALS S^NODFlT ■ TRUE* 
FINIS » FALSE* 
♦♦♦♦♦PLOT PARAMETERS****^ 

NX • 60* NY « <.5* 2MAP • 0* 1* 2* 3* 4* 5* b, 7* 8, 9, 

XMIN • 0, XMAX s 0* YHIN " 0* TMAX ■ 0, 

ZMIN • 0* 2MAX » 9, XRICH » 0* YKiCrt ■ 0* 

DEFAULT » 0* 0VPRfi!7 « F* AVE ■ F* I.^TZJ ■ ft 



/t**^^**'¥*^******'*^***if**i>.[jH l^f *♦«♦♦♦♦♦♦«**«***♦««* 

A - 5.<»8» & • 0.0'v* C ■ :?.67* S 

TAiJLE l.A CDNVERGtNCE OF THE ITERATED PARAMfcTEt?S 



ITERATITM 



10 
20 
30 
AO 
50 



<.,379 

<». 973 
5.35.% 
5.676 



«0A3 
.0^5 
.0<.5 
.O^.'i 
.043 



6.27:; 
6.300 
6.146 
5.942 
5 .733 



RESIDUAL 
.O29E+02 
.52yE+02 
.b21c+02 
.525E+02 
.532E+02 



TABLE l.tJ DATA(A*YJ AND APPROXIMATION (F) 



X( J ) 



Y{ J) 



F( J) 



.300f:-''-0i 
.53UE+01 
.960E+01 
.lfaOE+02 
.3C0E+02 



.390E+02 
. 550E+02 
.575E+02 
.600fc+02 
.570E+02 



.436E*32 
.500E*02 
.p78t*-02 
.623E+02 
.553E+02 
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NONLINEAR LfcAST SJUARES FIT OF 

Y-C*(X+A)*EXP(-B*X)> (2) 
TO CQHO SALMON GROWTH JATA (1) 
A* t>.68 a* .0<> Ce ^.73 



.300 
X — 

6.226" 

6 
R 
0 

W 5.962 
T 
H 



A 5.6V7 
T 

E 



G 
/ 
G 
/ 

0 i.l69 

A 

Y 

) 

4. 90<i 



.712 

-X — ■ 



1.169 
■ — X — 



<i. 6A0 



<i. i76 



<i.lll 



3.900Y1 

X— 
.300 



22 



1.627 2.065 
— .--.-X— --<- — --X— ■ 

2222222222 
. 222 222 
222 222 



2.5^.2 
X— 



3.000 
X 



22 
2 



222 
22 



22 



22 



22 
22 



22' 



2 1 
22 



2 
2 



2 
2 



— X — — — — — — X — — — — — — X X 

.712 1.169 1.627 2.065 2.5^2 

OISSULVED OXYGEN (M&/L) 



3.D00 



SCALE FACTORS ■ X-AXISJ £+01 Y-AXISJ E+01 Z-AXISj f*00 
ZO-Z^- 0(-9)/ 1.000( 5)* 2.000(6A)* 3.000( 0}/ A.000( 0) 

25-29- 5.000( 0)» 6.0001 0)f 7.000( 0)* B.000( 0), 9.000( 0) 
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PROGRAM -INT6- READY FOR INPUT 

/ / 

/ ■ / 
OPTION " Z, A ■ 0.70> B • 51. 0# i 

TABLE 2. A CONVERGENCE OF ITERATIONS OF L.S. INTEGRAL 

1220.205 
1213.685 
1213.665 
1213.665 
1213.665 
1213.665 



TAdLE 2.B 

VALUES FOR DATA (X/YJ* THE EXPONENTIAL FJNCTlUN 6lX)i 
AND THE LINEAR APPROXIMATION L(X) 

X Y G(X) L(X) 

3.00 39.00 <.2.6'r ' 53.10 

5.30 55.00 <.9.50 5<..71 

9.60 57.50 57.82 57.72 

16.00 60.00 62. Al 63.60 

30.00 57.00 55.35 72.00 
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LEAST SQUAiiES FIT OF A STRAIGHT LINE* (3)# 
TO Y-b.i'CX+S.O^EXPC-O.O^S^X), (2) 
WITH OP.IGINAL OATA> INCLUDED. 
LINE IS Y« .70*X ♦ 51.00 



.300 
X— 



6 
R 
0 
W 
T 
H 

R 
A 
T 
E 

( 
M 
G 
/ 
G 
/ 
0 
A 
Y 
) 



7.158 



6.<.18 



6.0^7 



5.677 



5.307 



A. 937 



^.566 



A. 19b 



3.900Y1 

X — 
.300 



► 712 
-X— 



l.lb9 
. — X— 



1.627 
• — X — 



2.085 
■ — X — 



2.542 
■ — X — 



33 



333 
333 



33 



3.000 
— X 
33 



33 



13i3 
333 
33 



3 33 
333 

33 

333 
333 

33 

333 
333222 
22233 22222 
222333 2222 
22 33 22 
22 333 1 222 

333 22 
233 22 
333 22 
3332 22,^1 
3322 22 
2 
2 



X X X--- 

.712 1.169 1.627 2.0tf5 

OISS'JLVEO OXYGEN (MG/L) 



• — X- — 
2.542 



• — X 
3.000 



SCALE FACTORS - X-AXIS* E+01 
Z0-Z4- 0(-9), l.OOOt 4)» 

25-29- 5.000( 0), 6.000{ 0), 



Y-AXISJ E+01 
2.000(54)/ 
7.000( 0)p 



Z-AXIS> 
3.000(67)* 
d.000( Q)f 



E+00 
4.000( 
9.000 I 



0) 
0) 
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PROGRAM -INT(i- READY FOR INPUT ^ 
/ / 

/ • ' 
OPTION • 3* ANGLE^" 1.0» ^OIAM • 1.0* 

.SOLAR ENERGY OIlStI ON^OUE '-TO ANGLE OF INCLINATION 

INCLINATION ANGLE ■ 1.00 HOLE DIAMETER • 1.0 

AREA OF CIRCULAR HOLE • .79 

AREA ILLUMINATED • 1.^5 DILUTION CUEFFICIEnT • .5A03 
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SDLAR OILUTIUN 
ELLIPTICAL REGION ON 6KJUN0 ILLUMlNATtO 
BY SUNLloHT PASSING THROUGH A CIRCLE OF 
OlAMtTER " 1.00, INCLINATION ANGLE • 1.00 RADIANS 



1. 



-1.851 

X- 

OOOY 
I 



-l.2bb 
X— 





X 


111 


0 


I 


11 


s • 


I 


11 


I 


• 773y 


11 


T 


I 


1 


I 


I 


11 


□ 


I 


1 


N 


I 


1 




• 5<i!> Y 


1 


\ 


1 


1 


L 


I 


1 


□ 


I 


1 


H 


I 


1 




• 3 1 d Y 


1 




T 

Si 




K 


I 


1 


I 


1 


1 


N 


I 


1 


0 


.091Y 


i 


R 


11 






11 




A 


U 




X 


I 


1 


I 


-.136Y 


1 


s 


• I 


1 




I 


1 


0 


I 




f 


I 


1 




-.3b^Y 


1 


E 


I 


1 


L 


I 


1 


L 


I 


1 


I 


I 


1 


P 


-.t>91Y 


I 


S 


I 


1 


r" 

c 


I 


11 



-.659 -.031 
._x — — — X — — — X — — ' 

111111111111 
11111 11111 

111 



1.223 
X — 



1.851 
X 



li 



11 



11 



1 
1 



1 
1 



1 
1 



1 
1 



1 

1 



1 
1 
1 



1 
1 



I 
I 

-. 3iaY 
I 
I 

I 

-l.OOOY 



-1.851 



11 

1 1 

\i 

11 11 
Hi 111 

11111 illll 
111111111111 

_-X X X X X — 

1.2ci6 -.659 -.031 .596 1 .223 

POSITION ALONG MAJOR AXIS OF ELLIPSE 



— X 
1.351 



SCALE FACTORS ■ X-AXiSl E+00 Y-AXISi fc+OO Z-AXISi E ^00 
fO-Z'." 0(-9)# 1.000(-9), 2.000( 0)* 3.000( 0)^ 't.OOCl 0) 

i'j-Z=)' 5.000( 0)» 6.000( 0), 7.000( 0), 8.000( 0)* 'V.OOCl 0) 
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PROGRAM -INTG- READT FUR INPUT 

'> ; 

♦»4,»**4SYQp pROGRArt**»*»»»*»»»»»»»*»»*»*»»**/ 

FINIS ■ TRUE#S ' 

PROGRAM -INTG- TERMINATED 
EOR 
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